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Abstract 

We discuss an effective Polyakov loop model for QCD thermodynamics 
with a chemical potential. Using high temperature expansion techniques the 
partition sum is mapped exactly onto the partition sum of a flux model. In the 
flux representation the complex action problem is resolved and a simulation 
with worm-type algorithms becomes possible also at finite chemical potential. 
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Introductory remarks 



With the running and upcoming experiments at BNL, CERN, GSI and Dubna the 
amount of experimental facts about the QCD phase diagram will increase consider- 
ably in the near future. Also the theoretical side is challenged to contribute to our 
understanding of QCD thermodynamics, a task which is rather demanding due to 
the non-perturbative nature of the problem. In principle lattice QCD is a suitable 
non-perturbative approach, but at finite chemical potential the complex phase prob- 
lem considerably limits the applicability of Monte Carlo methods for a numerical 
evaluation of the path integral. 

A tempting idea to overcome the complex phase problem is to search for a 
transformation to new degrees of freedom where the partition function becomes a 
sum over configurations with real and positive weights. Such representations with 
real and positive weights can then be used in a Monte Carlo calculation. Several 
examples for this kind of transformations in various models can be found in the 
literature, but for QCD we are still far from identifying suitable degrees of freedom. 

In this article we contribute to the enterprise of developing new and real repre- 
sentations for QCD related systems with a quark chemical potential fj, by considering 
an effective theory for Polyakov loop variables. The Polyakov loop has the inter- 
pretation of a static color source and in pure gauge theory its expectation value 
serves as an order parameter for the deconfinement transition. More specifically 
it is an order parameter for center symmetry which is broken spontaneously at the 
deconfinement transition of pure gauge theory. Although in full QCD the underlying 
center symmetry is broken explicitly by the quarks, also there the Polyakov loop may 
be used to determine the crossover to deconfinement. 

The idea of using effective theories for the Polyakov loop to describe the transi- 
tions in pure gluodynamics and in QCD goes back to the late seventies, and since 
then, both the understanding of sub-leading terms in the effective action, as well as 
the simulation techniques have been improved considerably [1]-[10]. 

The effective theory for the Polyakov loop considered here may be derived from 
QCD using strong coupling expansion for the gauge action and hopping expansion 
for the fermion determinant. The leading center symmetric and center symmetry 
breaking terms are taken into account, as well as the chemical potential /i. For 
fi the action becomes complex, i.e., in its standard representation the theory 
inherits the complex phase problem from QCD. 

In this work we apply high temperature expansion techniques and derive a flux 
representation which is free of the complex phase problem. The new degrees of 
freedom are integer valued flux variables attached to the links of the lattice and 
integer valued monomer variables at the sites. The flux-monomer configurations are 
subject to constraints which enforce the total flux at each site to be a multiple of 
three. The weight factors for the admissible configurations are given in closed form 
and turn out to be real and positive. Thus, in the flux representation a Monte Carlo 
simulation of the system with generalized Prokof 'ev-Svistunov worm algorithms [llj 
becomes possible and should allow for an effective numerical treatment without any 
complex phase problems. Finally, we show how observables can be expressed in 
terms of the flux and monomer variables. 
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The effective theory and its high temperature expansion 



The effective Polyakov loop theory we consider is described by the action 

3 

S[P] = ~t^2Y^ [TrP(a;)TrP(a; + z>) t + TrP(» t TrP(x+i>) 

X U — l 

- [??TrP(x) + rjTcP(x)^ . (1) 

X 

In this standard representation the degrees of freedom are the SU(3) valued Polyakov 
variables P(x) attached to the sites a; of a three-dimensional cubic lattice which 
we consider to be finite with periodic boundary conditions. By v we denote the 
unit vector in ^-direction, with v = 1,2,3. The first term of the action, which may 
be obtained from a strong coupling expansion, is a nearest neighbor interaction of 
the traced Polyakov loops. This expansion also establishes r to be an increasing 
function of the QCD temperature, and for simplicity we refer to r as temperature. 
The term in the second line of ([TJ can be derived from the fermion determinant 
using hopping expansion. The parameters rj and rj are related to the amplitude k 
and the chemical potential (i via 

X] = ne^ , rj = Ke~^ . (2) 

The hopping expansion identifies n to be proportional to a power of the hopping 
parameter and thus n is a function of the quark mass m which decreases with 
increasing m. For Nf flavors of mass degenerate quarks k is simply proportional to 
Nf. The terms e^TrP(x), e~ fl TrP(x)^ in {TJ, fl2J are the leading /it-dependent 
terms of the fermion determinant 

From Eq. (TJ} one immediately sees that the effective theory has a complex 
phase problem: For /i ^ one has r) ^ rj, which leads to a non-vanishing imaginary 
contribution to the action given by i^2 x (fj — rj)lmP(x). Thus the Boltzmann 
factor exp(— S[P}) obtains a complex phase and cannot be used as a probability 
weight in a Monte Carlo calculation. 

Effective theories of the type ([TJ were studied before in various contexts. For 
the case of vanishing n an important line of research is the determination of sub- 
leading terms in the action to get the effective theory ready for quantitative analysis 
- see, e.g., [Z] _ [8] For the case of non-zero chemical potential an important 
contribution is [10J, where the theory described by ([TJ was studied with complex 
Langevin techniques. A future numerical simulation with the flux representation 
presented here should be able to check the results [10J and shed light on the reliability 
of the complex Langevin method. 

A reduced version of the model (pj, where the traced Polyakov loops are re- 
placed by center elements, was studied as well |3J-|6J. In this simpler case a flux 
representation free of complex phases has been known for a long time |3J and various 
numerical simulations with different techniques were presented |3J-[6J. 



1 Actually in the fermion determinant the chemical potential appears rescaled with the 
temporal extent Nt of the lattice, but in order to simplify the notation used for the effective 
theory we omit this factor. 
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The grand canonical partition function of the model described by JTJ is obtained 
by integrating the Boltzmann factor e~ s ^ over all configurations of the Polyakov 
loop variables. The corresponding measure is a product over the reduced Haar 
measures dP(x) at the sites x. Thus 



Z = Jl[dP(x) e- s ^ = jD[P}e~ s ^ 



(3) 



where in the second step we have introduced the shorthand notation D[P] = 
Y[ x dP(x) for the product measure. 

The first steps towards the flux representation are writing the Boltzmann factors 
in a factorized form and an expansion of the remaining exponentials. This step 
corresponds to an expansion in r. In spin system language r would have to be 
identified with the inverse temperature f3 and thus our expansion is equivalent to 
high temperature expansion in statistical mechanics. We obtain the following form 
of the partition sum: 



Z= fD[P][l[e T 



Trf(a;)TrP(i+P) f £ t Tr P(k) 1 " Tr P(x+0) 




^TrP(x) WTrP(x) + 
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(4) 



In the end the partition function will be a sum over configurations of the expansion 
coefficients £ a ,„,Z X)I , for the nearest neighbor terms, and the expansion coefficients 
s x ,~s x for the site terms. We will refer to l XM: l x , v as flux variables and to s x ,s x 
as monomer variables. For the sum over flux- and monomer configurations we 
introduce the shorthand notation 




(5) 



Using this notation, rearranging the products in (j4|), and writing the integral over all 
configurations of the Polyakov loops P(x) again in its factorized form, the partition 
sum reads 



*= e n 



{1,1, s,s} 



~t\ \XM 



I I / I 



n 



T) Sx T) Sx 



(6) 
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J| /dP(x) (TrP(a;)J^ w (TrP^) 



x.v ~T" lx — i>,v\ S x \ / j [lx,u l>x — u,u\ ~t~ 

In the last step we introduced the weight factor for a configuration of flux- and 
monomer variables, 

- (g SS) (? 55) ■ (7) 

and for the remaining SU(3) integrals at the sites use the abbreviation 

I(n\m) = / dP (TtpY (TtpA™ . (8) 



Here n and m are non-negative integers and JdP denotes the integration over 
SU(3) Haar measure. 



Solving the SU(3) integrals 

In this section we evaluate the remaining SU(3) integrals I(n\m). The corresponding 
generating function is the integral 



G(u\v) = dPe uT * p e vTlP ' = £ I{n\m) , (9) 



iTrP „i)TrP t _ 

n,m— 

in other words the I(n\m) are the moments of the one link integral G(u\v). We 
start our derivation of the I(n\m) with an expression for G(u\v) given in |12j , 

G{u\v) = E r + +lVr 2 ++2V , ( 3(P+g + 1) ) («»)'(«-+«»)«. 

We use the binomial formula to evaluate (u 3 + v 3 ) q and organize the terms with 
respect to the monomials u 11 v m to obtain, 

? f 3 (P + 9 + l)^ \^ ( <l\ u P+3j v P+3 q -3j 



g(u\v) = y - ^ y , 

n! m! ^ (p+ q + 1)\ ( p+q +2)1 q\ 

n,m=0 p,q=0 j=0 ^ H ) \f ' H ' / 1 

Comparing this expression with the second form of the generating function ([9]), we 
identify 

~ « 2n\m\ (3(p+q+l)\(q\ 
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The two Kronecker deltas can now be used to reduce this expression to a finite sum. 
The first one implies p = n — 3j, and since p > 0, we find n > 3j and obtain an 
upper bound for the sum over j given by j < [n/3\, where by [x\ we denote the 
floor function^]. Thus after performing the sum over p we have 

2n\m\ f 3(n— 3j + q+l)\f q 

~ tho f^> ( n _ 3j+q+1 y( n _ 3j+q+2)lql y n _ 3j 

The remaining Kronecker delta may be written as <5n- m ,6j-3<?. which makes explicit 
an important property of the integrals I(n\m), the triality constraint: 

I{n\m) ^ , only if (n — m)mod3 = 0. (10) 
The Kronecker delta Sn-mfij-Sq implies q = 2j — (n — m)/3. Using this in the 
second binomial factor gives = ^ 2 J~("7 m '/ 3 ^ . In order to obtain a non-zero 
result for this binomial, the upper argument may not be smaller than the lower one, 
i.e., 2j — (n — m)/3 > j. This implies j > (n — m)/3, and since j also must be 
non-negative, we obtain the lower bound j > max(0, (n — m)/3). Performing the 
sum over q and using the lower bound for j we obtain the I(n\m) as finite sums 

I(n\m) = V - - ^ — ^ - '— 



j— max (0, 



Here we have introduced the triality function 



3 

(11) 



, . J 1 for nmod3 = 0, , 
T M = { else . < 12 > 

The final result (fTT|) expresses the moments I(n\m) as finite sums. In [12J recursion 
relations that relate different I(n\m) were presented and in an appendix the lowest 
moments for n,m < 10 are listed. We compared the results from our explicit 
expression (jlip to these values from the recursion relation and found agreement. 



Final form of the partition sum and graphical representation 

Using the results from the last section we obtain our final expression for the partition 
sum in terms of flux and monomer variables, 

Z= ^ C[U,s,s] • (13) 

{1,1, s,s] 

The first term under the sum is the total weight factor assigned to a configuration 
of fluxes and monomers, 

mJ,s,s] = [n^^^fnTrlr)^ 1 ' 5 ^] with 

\ „. „ fx, i/! lx,v- / \ _ b x- b x- ) 



' \ x\ is the integer with x — 1 < \x\ < x. 
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X[lJ,S,s] = ]Jll^[l x ,u+L-u,u] + S x ^2\l x ,u + lx-0,u} + S x \. (14) 

The second factor under the sum in ([TBI is the constraint, 



C[1,1,S,S] = n T (l] [(lx,u ~lx,u) ~ (J X 



(Sx 



(15) 



The constraint is a product over the sites x and at each site the triality function T 
enforces the combination 



fx 



E 



^x.u) {J-x — um ^x — Om) H~ {^x $x) 



(16) 



of flux and monomer variables to be a multiple of 3. The expression f x is the total 
net flux at site x, including the contributions from the monomers. 

It is instructive to compare the flux representation (|13p - (|15p to the flux rep- 
resentation [2, 6J for a simpler effective theory where the Polyakov loop variables 
TrP(a;) are replaced by Z 3 center elements P(x) G {1, e l2jr/3 , e~ l2Tr / 3 }. The main 
difference is that in the discrete Z3-case for the dimer and monomer variables only 
the three values —1,0,-1-1 are necessary. Here, where the dynamical variables are 
TrP(x), i.e., continuous degrees of freedom, integers with an infinite range are 
necessary for the monomers and dimers. The triality constraint for the flux f x at a 
site x is essentially the same in the two models, with the difference, that for the Z3 
case again only a finite number of values is possible, f x G {—6, —3, 0, 3, 6}. Finally, 
for the effective theory studied here, the weight (Tl4|) contains the factors I(n\m) 
which reflect the non-abelian nature of the degrees of freedom TrP(x). 

We now present a graphical representation for the flux and monomer variables 
an discuss the triality constraint for the net flux f x in this graphical language. The 
flux variables l xv and l x v are assigned to the links of the lattice. In Fig.[TJwe show 
a site x, as well as its two neighbor sites along the v direction, and indicate the 
flux variables assigned to the corresponding links. The monomer variables s x and 
s x are located at the sites, and again Fig. [TJ illustrates the assignment. 



^X-V f &X-V , Sx f Sx 



Sv+V * S 



X+VJ "X+V 



lx-v,v > lx-%v I lx,v > lx,v i 



x-V 



X 



X+V 



Figure 1: Assignment of the flux variables l Xt v,l x ,v, and the monomer variables 
s x , Sx to the links and sites of the lattice. For clarity we display only one of the 
three possible directions. 
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Figure 2: Examples of admissible flux- monomer configurations at a single site 
x. For clarity we show only configurations in the 1-2 plane (see the axis labels). 
A flux variable of l x>v = n is represented by n arrows pointing in positive in- 
direction, while the l x>v point in negative ^-direction. A monomer variable 
s x = n is represented by n outwards pointing (i.e., pointing away from the site 
x) triangles, while for the s x we use inwards pointing triangles. The details of 
the four examples a) - d) are discussed in the text. 

In our graphical language we use arrows and triangles to represent the flux and 
monomers (see Fig. [2] below). For a flux variable l XtV = n we put n arrows pointing 
in positive v direction on the link between x and x + v. For l xu = n we use n 
arrows pointing in negative ^-direction, since the l x>v enter with a relative minus 
sign in the expression (fl6|) for the local flux f x . Similarly, for a monomer value 
s x = n we use n outwards (i.e., away from x) pointing triangles, while the ~s x 
are represented by inwards pointing triangles. We stress that Fig. [2] shows only 
flux-monomer configurations for a two-dimensional sub-lattice, the 1-2 plane, and 
the vertical dimension in the plot is used to display the triangle symbols for the s x 
(pointing outwards) and s x (pointing inwards) monomers. 

In Fig. [2]we show four examples of admissible flux-monomer net fluxes f x at a 
site x. The first example (Fig. [2^) is a double line of flux entering from the negative 
1-direction, which then splits at x and two single fluxes exit in 1- and in 2-direction. 
The non-vanishing variables are l x _i 1 = 2,/ Xi i = 1,^,2 = 1- All other flux and 
monomer variables attached to x vanish. The total net flux is f x = 0. In Fig. [2p 
a double line of Z-flux enters and is compensated by two units of the monomer 
variable {l x _2 2 = 2> s x = 2, f x — 0). Fig. [2]: shows an example of a non-zero 
total flux f x = —6 generated by a combination of nonvanishing I- and 7-fluxes, as 
well as monomers {l x _i A = 1,^_2. 2 = = Mx,i = 2,1^2 = l,s x = 2). 
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Finally Fig. [2jd is an example with f x = —3, described by I j x — 3, l x _§2 = 
2,'x,i = 1)^,2 = = 2,I X) 2 = l,Sa; = 1,^ = 2. The configurations in the 

sum J2{ilss} are * nen obtained by combining admissible flux arrangements at all 
sites x. 

The graphical representation will be useful for developing generalized Prokof'ev- 
Svistunov worm algorithms [llj which may be used to update the flux representation 
efficiently |13| . Furthermore, one can expand the partition sum (fl3|) for small r. 
Such a result will be important to check the outcome of numerical simulations (see 
also |6J), and the graphical representation is a highly welcome tool for organizing 
the terms of such an expansion. 



Flux representation of observables 



For a successful application of the flux representation of the effective theory in a 
Monte Carlo simulation also the observables have to be expressed in terms of flux and 
monomer variables. We here briefly discuss this issue for simple bulk and fluctuation 
observables. In particular we focus on the expectation value (P) of the Polyakov 
loop, where P = J2 X TrP(x), the corresponding susceptibility \p = {P 2 ) ~ (-P) 2 . 
the internal energy U = (S) (where S is the action {!])), and the heat capacity 
C = (S 2 ) — U 2 . These quantities may be obtained as simple derivatives of the 
partition sum, 



(17) 



Xp 



U 



Z 
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z 



° = -z 
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V 
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drdi] 



+ 2rr/ 



drdfj 



d7]dfj 



Z-U 2 



From the final form JT3J) of the partition sum it is obvious that the derivatives in 
(fT7|) may be pulled under the sum over all flux and monomer configurations and 
there only act on the weight factor W[i,I, s,s], since the constraint is independent 
of r, rj and rj. To perform the needed derivatives we write the weight factor as 



yv[i~i,s,s] = t l+l v s t] s (\ 




l[l,l,s,s], (18) 



where we introduced abbreviations for the sums of flux and monomer variables 

L = ^^l x ,u , L = *^2~l x ,v , S = , S = ^s x . (19) 
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Using the weight factor (|18p one immediately obtaines the following relations: 



£-W[l,l,8,?] = -W[l~l,s,s] 
or] rj 



s 2 -s 

rf 



W[l~l, 



s,s\ 



(20) 

for the derivatives with respect to rj, and similar relations for the derivatives with 
respect to the other parameters. The important aspect of the expressions (|20|) is 
that the derivatives of the weight factor may be written as the unchanged weight 
W[l,l,s,s], multiplied with factors built from the parameters 77,77, r and the sums 
(fl9|) of flux- and monomer variables. Reinserting these relations we obtain rather 
simple expressions for our observables (jTTJ) , 



1 



(21) 



XP = 



U = (L-\ 
C = ( ( L 



(S 2 -S) - (S) 1 
L + S + S) , 

f I+S + S) 2 ) - U - U 2 



Obviously all the bulk and fluctuation observables we list here can be expressed in 
terms of the (summed) flux and monomer occupation numbers from Eq. (|19f> , as 
well as the corresponding fluctuations. 

As a matter of fact, in a similar way also correlators can be expressed in terms 
of the flux and monomer variables. We briefly discuss this for the example of the 
Polyakov loop correlator ( TrP(x) TrP(y) t ). The two factors TrP(x) and TrP(y) t 
can be generated from the partition sum by using parameters rj(x) and rj(x) in ([!}, 
which are independent parameters for all sites x of the lattice. The mapping of 
this generalized model with locally varying parameters to the flux-monomer repre- 
sentation goes through unchanged. The local parameters 77(2) and rj(x) can now 
be used as sources and the correlator ( TrP(a;) TrP(y)^ ) is obtained via two partial 
derivatives. In the end one sets all parameters to the original values, i.e., rj(x) — r\ 
and rj(x) = rj for all x. One ends up with 



(TrP(cc)TrP(y) t ) 



O 2 



Z drj{x) drj(y) 



(22) 



We find that the Polyakov loop correlator turns into a correlator of the monomer 
numbers at the corresponding sites. In a simulation of the flux-monomer represen- 
tation with a worm algorithm one can evaluate this correlator by directly sampling 
the starting point of the worm correlated with its actual position and with this 
improved estimator one expects to get an excellent signal. 



Summary, discussion, outlook 

In this article we study an effective Polyakov loop theory for QCD thermodynamics. 
The model may be obtained from a strong coupling expansion, combined with 
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a hopping expansion for the fermion determinant. It contains the leading center 
symmetric and center symmetry breaking terms. At non-zero chemical potential the 
model inherits the complex phase problem of finite density QCD. The complex phase 
problem prohibits a direct Monte Carlo simulation in the standard representation. 

Using high temperature expansion techniques, we map the partition sum of the 
effective Polyakov loop model to a representation where the dynamical degrees of 
freedom are integer valued flux variables attached to the links of the lattice, and 
integer valued monomer variables at the sites. The flux-monomer configurations are 
subject to a constraint which forces the local net flux at each site to have vanishing 
triality, i.e., the flux has to be a multiple of three. All admissible configurations come 
with a real and positive weight factor, and thus in the flux-monomer representation 
the complex phase problem is solved. 

For systems with such a flux representation generalizations of the Prokof'ev- 
Svistunov worm algorithm [llj allow for a numerical Monte Carlo simulation of the 
model. As a matter of fact, for a similar QCD-related system a generalized worm 
algorithm was used to effectively explore the full temperature and chemical potential 
parameter space, also in regions where the complex phase problem is severe |6J. 
We expect [1 3J that for the system discussed here, a worm algorithm with similar 
efficiency can be implemented on the basis of the new flux-monomer representation 
given here. Again it should be possible to explore the full temperature-density phase 
diagram. Such an analysis will provide interesting insight into various aspects of 
the phase structure of QCD-related systems. 

Besides the improvement of our physical understanding, a simulation based on 
the new flux representation will also have important technical applications. The 
models can serve as a prototype system, which other approaches to QCD with non- 
zero chemical potential can be tested against. Various expansions in fi, imaginary 
chemical potential techniques, as well as results from reweighting (see [14] for 
reviews) could be compared to the outcome of a Monte Carlo simulation in the flux 
representation, where the simulation is not tainted by uncontrolled effects from a 
complex phase problem. 

However, not only lattice techniques should be compared to the prototype system 
presented here. It is straightforward to write down a continuum counterpart of 
the action JTJ), and the system could also be analyzed with various continuum 
techniques, in particular functional methods [15j. A comparison of these results with 
the lattice would provide important insights on the assumptions and the reliability 
of the continuum techniques. 
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